# 渗流模型模拟程序

import numpy as np
import sys
# 设置递归深度为100000
# 因为本程序多次使用for循环和嵌套递归，默认情况下递归深度达到1000会出现异常
sys.setrecursionlimit(100000) #这里设置为十万

# 程序分为四大模块
# 一、矩阵生成模块
# 二、初始值定义模块
# 三、团簇计算模块（主函数模块）
# 四、结果输出模块


# 一、矩阵生成模块
# 矩阵生成函数
def new_array(array,rate):  #array是定义好的全1矩阵，rate是矩阵中每个值取0的概率
    '''按照数组模板生成对应的 0-1 矩阵'''
    zeros_num = int(array.size * rate)  #根据0的比率来得到 0的个数
    new_array = np.ones(array.size)  #生成与原来模板相同的矩阵，全为1
    new_array[:zeros_num] = 0  #将一部分换为0
    np.random.shuffle(new_array)  #将0和1的顺序打乱
    re_array = new_array.reshape(array.shape)  #重新定义矩阵的维度，与模板相同
    return re_array

# 生成初始0，1矩阵
P = float(input('请输入概率P：'))
size = int(input('请输入矩阵的长（宽）：'))
arr = np.ones(size*size)  # 生成一维数组arr，值默认为1
arr2 = arr.reshape(size, size)  # 将数组array变形成15*15的矩阵arr2
new_arr = new_array(arr2,rate=P)  # 设定每个数字变为0的概率 rate


# 输出生成的初始0，1矩阵
print('----------初始0，1矩阵为----------')
print(new_arr)
print('--------------------------------')


# 生成全0 N矩阵
# 该矩阵是为了统计在几行几列的位置上存在最大、次大团簇
N = np.zeros((size , size))


# 设置最大团簇计算矩阵（N）索引的初值
l = 0  # l是上面生成的矩阵N行数的索引控制值
k = -1  # k是上面生成的矩阵N列数的索引控制值
M = []  # M为了容纳之后矩阵N的列表转化

def up(i,j): # 负责控制团簇向上方向遍历
    global l
    global k
    if i == -size - 1: # 因为矩阵的遍历范围只有（-size,size-1）为了防止遍历时超出索引范围，需要重置索引值
        i = size - 1
    if new_arr[i][j] == 0:  # 当遍历到的值为0时，说明该点为团簇
        N[l][k] = N[l][k] + 1  # 每发现一个团簇便记录一次团簇大小
        new_arr[i][j] = 2  # 每当记录完一个团簇便重置该点的值为2，以防之后的重复计算
        e = i
        f = j - 1
        left(e, f)
        g = i
        h = j + 1
        right(g, h)
        i = i - 1
        up(i, j)
    else:
        N[l][k] = N[l][k]


def down(i,j): # 负责控制团簇向下方向遍历
    global k
    global l
    if i == size:
        i = -size
    if new_arr[i][j] == 0:
        N[l][k] = N[l][k] + 1
        new_arr[i][j] = 2
        e = i
        f = j - 1
        left(e, f)
        g = i
        h = j + 1
        right(g, h)
        i = i + 1
        down(i, j)
    else:
        N[l][k] = N[l][k]


def left(i,j): # 负责控制团簇向左方向遍历
    global k
    global l
    if j == -size - 1:
        j = size - 1
    if new_arr[i][j] == 0:
        N[l][k] = N[l][k] + 1
        new_arr[i][j] = 2
        a = i - 1
        b = j
        up(a, b)
        c = i + 1
        d = j
        down(c, d)
        j = j - 1
        left(i, j)
    else:
        N[l][k] = N[l][k]


def right(i,j): # 负责控制团簇向右方向遍历
    global k
    global l
    if j == size:
        j = -size
    if new_arr[i][j] == 0:
        N[l][k] = N[l][k] + 1
        new_arr[i][j] = 2
        a = i - 1
        b = j
        up(a, b)
        c = i + 1
        d = j
        down(c, d)
        j = j + 1
        right(i, j)
    else:
        N[l][k] = N[l][k]


def Sagiri(): # 团簇遍历的总控制函数，负责综合矩阵中四个方向的遍历，确保矩阵中的每个值都得到遍历
    global l,k
    for i in range(0, size, 1):
        for j in range(0, size, 1):
            l = i
            k = j
            if new_arr[i][j] == 0:
                N[l][k] = N[l][k] + 1
                new_arr[i][j] = 2
                a = i - 1
                b = j
                up(a, b)
                c = i + 1
                d = j
                down(c, d)
                e = i
                f = j - 1
                left(e, f)
                g = i
                h = j + 1
                right(g, h)

            if new_arr[i][j] == 1:  # 当矩阵中值为1时，不计入团簇，并结束该次循坏
                continue
Sagiri()


# 四、结果输出模块
# 输出团簇分布情况的矩阵N
print('----------团簇大小分布矩阵为----------')
print(N)
print('-----------------------------------')

# 将矩阵N转化为列表N，用于对团簇大小进行排序
for i in range(0, size, 1):
    for j in range(0, size, 1):
        M.append(N[i][j])
A = (np.sort(M))


# 输出团簇大小排序
print('----------团簇从大到小排序为----------')
print(A[::-1])
print('-----------------------------------')


# 计算团簇的平均大小
Sum = 0
Ave = 0
s = 0
for i in range(size*size-1,0,-1):
    if A[i] != 0:
        Sum = Sum + A[i]
        s = s + 1
        Ave = (Sum)/s
    else:
        break

#输出团簇的平均大小
print('----------团簇的平均大小为------------')
print(Ave)
print('-----------------------------------')